From microfacets to participating media: A unified theory of light 
transport with stochastic geometry — Supplemental 


DARIO SEYB, Dartmouth College, USA 
EUGENE D’EON, NVIDIA, New Zealand 
BENEDIKT BITTERLI, NVIDIA, USA 
WOJCIECH JAROSZ, Dartmouth College, USA 


1 DERIVING RECURSIVE ENSEMBLE AVERAGE LIGHT 
TRANSPORT 


The step from Eq. (12) to Eq. (14) in the main document is not 
immediately obvious. It does not require any deep insights, but 
a thorough reproduction of the intermediate steps, while lengthy, 
provides extra intuition for our method. We provide this here in a 
slightly more general form that should apply to any quantity that 
can be estimated with a recursive Monte Carlo estimator. 


1.1 Expectations over realizations of a Gaussian Process 


We will often want to compute expectations “over realizations of a 
Gaussian Process”, that is, functional integrals of the form 


Peun [ vay LE Sak (a) 


where f is a realization of the Gaussian process GP(p, k), £ an 
operator acting on f, and y,,,(f) the classical Wiener measure of 
f with respect to GP(, k), i.e. the probability density of sampling 
f ~ GP(p, k). If the operator £ is linear in f, we can resolve this in 
the following straightforward fashion 


f, ET Lf dyp) =L f, m fay P) =Le (2) 


but many operators, in particular light transport, which we study in 
this paper, are not linear in f, and the above simplification does not 
hold. In the following, we will drop the mean j and covariance kernel 
k for notational convenience. A rigorous treatment of functional 
integrals of this form is out of the scope of this work. Instead, we 
can provide intuition for our results by restricting ourselves to 
Gaussian processes over a finite index set X with |X| = n, and 
discrete operators £”. In that case, the Gaussian process is simply a 
multivariate normal distribution in n dimensions, and we can write 


1 £”fay(f) = 
GPx 
fon fete fal a" Ohi fal) = 
R R 
f[ L"Fdy"(F), (3) 
R” 


where y” is n dimensional Gaussian measure over X with the same 
mean and covariance as the given Gaussian process. Informally, we 
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write 
f Lf dy(f) = lim 1 L"F dy” (F). (4) 
GP noo Jpn 


Of course, there is no immediate reason to believe that this limit is 
well defined in any meaningful sense. That said, this technique, also 
known as “time slicing” [Grosche and Steiner 1998], is commonly 
used to make sense of functional integrals, and for the purposes 
of our work and the domains we consider, the time slicing view is 
sufficient. 

In particular, it lets us treat Gaussian processes like we treat finite- 
dimensional distributions in many situations. For example, we can 
write 


f Lfay( = f LL fan fae) Ay Ufa. far) 
GP GP 
. f f LL fa, fac] dyac (fae | fa) dya(fa) 
GPa JGP acl fa 
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(5) 
where A is a subdomain of the domain of GP, A‘ its complement, fA 
and fac realizations of GP over A and A‘ respectively and [ få, fac] 
the “concatenation” of the individual realizations. The first equality 
simply decomposes f into f4 and f4c. Going to the second line, we 
apply the law of conditional expectations. That is, we can sample 
f “all at once” or first sample a part of f (f4) and then the rest 
(fac) conditioned on the part that we sampled first. The last line 
simply says that we don’t have to restrict the inner process to A‘ 
and can instead “resample” over the whole domain. For the part of 
the domain covered by the “conditioning region” A, we will simply 
get back the same values that we conditioned on. 


1.1.1. Conditioned Expectations. Additionally, we can write condi- 
tional expectations as usual: 


(LL fac, al) Gp ac 
ya(a) 


1.1.2 Expectations of delta functions. We can write an expectation 
containing a delta function over some subset of the domain A as 


(S( f(A) — a) £ for = i „FA -LFA 


(LF cpl f(A)=a = (6) 


: f 5(fa - a) f Lf ay(f |fadyafay P 
GPa GP| fa 
= yalayX(Lf era 


In the first step, we simply expand the ensemble average. Next, we 
apply the decomposition in Eq. (5) and pull out 5( f(A) — a) from the 
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inner integral. We can then apply the definition of the delta function 
and drop the outer integral, replacing the integration variable fA 
with a and evaluating the measure y,4(f4) at a. 


1.1.3. Expectations of indicator functions. We can write an expec- 
tation containing an indicator function over some subset of the 
domain A as 


(f(A) > 0) £fop = f, ICA >O LFA) 
: I(f(A) > 0) ‘i Lf ay(f | fa) dya (fa) 
GPA GPac| fa 


E : f Lf dy(f | fa) dya (fa) 
GP% 9 GP| fa 
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where GP% is the process restricted to all positive realizations over 
the index set A. Note that A can be a manifold in the domain of the 
GP (e.g. a 1D line segment in a 3D domain). 


1.2 Application to light transport 
With all of the puzzle pieces in place, we can apply them to Eq. (12) 


from the main paper, letting Lf := L; (x, œ), we are particularly 
interested in the inner integral 


= f, forf 
i If p(x) i nce (xz, m) M0, t) Lf (xe, +) dyz (f) dos dn dt, 


(Sf (xen) 0,0) LI (xer) Ye 


(9) 
and pull it out for more compact derivations, treating n and t as free 
variables. We first apply Eq. (7) 


(8(f (x1) = 0) - AF f (xe) — n) - ÑO, t) + Lf (x, 04))¢ 


f F (10) 
= Yx,|7 (0, n)(V(0, t) - L; (xz, Ot) ey £(x,)=0AV f (xp =n? 


and then Eq. (8) and arrive at 
f 
¥x,|¢(0,m)((L; (Xt, 2) fox) ert) ICAF (xp )=0AV f(x, )=n° (11) 


Since f(9,x,) is already conditioned on {Af (x+) = OAV f (xt) = n, we 
can add these to the inner condition without changing the expected 
value. 


Yalt OD (Xr ONAL AF ax) Pay ery 0D 


where čs = f (xr) =0A Vf(x;) = n. Plugging this back into Eq. (9) 
and expanding the outer ensemble average into its integral form, 
we get Eq. (14) from the main text. 


1.3. Deriving the GPIS density in the Renewal and 
Renewal+ models 


In Section 5 of the main document, we use the GPIS density to make 
connections to microfacet and participating media theory. We briefly 
want to give the full derivation that connects Eq. (14) and Eq. (17). 

Recall that for the Renewal and Renewal+ models, we only con- 
dition on values at path vertices, not values along path segments. 


That means that ¢ k JR does not contain f(o,x,) and hence does not 
depend on the integration variable in the innermost integral. This 
allows us to pull the recursive ensemble average out of the inner 
integral in Eq. (14) of the main text as 


(Li(x", @))z 
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(xxt) 


T(x: 1%) 


Replacing the now “empty” inner integral with T(x; | ¢) gives us 


(Li(x", w) 


=f IESE 0, n | HTx: | Li (xt, OL) DEAL e da; dn dt. 
(15) 


And finally, defining T (t,n | ¢) := yx,(0.n | O)T(x | 2), we get 
Eq. (17) in the main text. 


2 GAUSSIAN PROCESS DETAILS 
2.1 Kernel Functions 


We give analytic forms of a set of kernels and their spectral densities. 
Note that our method is not restricted to only these kernels since 
we support any smooth kernel. We pick these kernels in particular 
because they cover a range of different fundamental properties, such 
as non-locality and non-positiveness, which have the potential to 
strongly affect the results of our method. We give their most simple 
forms and denote the common parameters standard deviation and 
length scale with o and I, respectively. Deriving the spectral density, 
and especially sampling from it, is often non-trivial and depends 
on the number of dimensions. We give spectral densities for d = 3 
where we can, and d = 2 otherwise. While theoretically, one does 
not have to sample basis functions proportional to the spectral den- 
sity (one could choose a set using some proposal distribution and 
re-weight, like in any other Monte Carlo estimator), sampling pro- 
portional to the spectral density drastically reduces the number of 
basis functions required. For a visual overview of the listed kernels, 
see Fig. 3. 


Squared Exponential. This is probably the most widely used kernel 
in Gaussian process literature. 


lx -= yll? 


k(x, y) = o exp- z 


(16) 
oe ilo? 


[a 
E 


We use this as the “default” kernel for our method. It produces 
smooth realizations and computations with it are very numerically 


p(w) = (17) 
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stable. It also induces the most “exponential” first-passage times 
and is a great choice when the aim is to represent close-to-classical 
participating media. 


ALGORITHM 1: Sample p** 


ë ~ U(0,1) 
r — y2 - log(é) 
h ~ U(0, 27) 


return [sin(ġ), cos(ġ)]" -r 


Rational Quadratic. The rational quadratic kernel is a superposi- 
tion of infinitely many squared exponential kernels with different 
length scales. The weighting of these length scales is determined by 
the additional parameter a. 


R lx- yl y? 
EL, y) = 02 h: $ ai (18) 
zási 
zito (4) * * lonk, (2) 
RQ Vaz 
Pa (@) = (19) 


In particular, limg soo ka (x, y) = k€ (x, y). For small a, the realiza- 
tions of the rational quadratic kernel have “fractal” properties. 


ALGORITHM 2: Sample p°2 
t~T(a,l-*) 


return Sample p*= with 1 = r7? 


Periodic. Periodic kernels produce periodic realizations. This strongly 
affects the light transport in the scene and is a very challenging case 
for our algorithm. Any limit to the memory will produce drastically 
different results. Additionally, this shows that the covariance matrix 
is not sparse in general. 


2 2 —y|lA) (20) 


(21) 


Locally Periodic. To ease some of the difficulties that come with 
the periodic kernel, while still preserving its interesting properties 
locally, we make use of kernel composition to construct a “locally 
periodic” kernel. 


Ke (x, y) = o° exp |- 


KLO Per (x. y) = K(x, y) . KE (x, y) 
pPoO) = (ph * p(w), 


where « denotes convolution. This kernel is “locally periodic” in the 
sense that in realizations new each “repetition” is allowed to deviate 
slightly from the previous one. The length scale of the squared 
exponential kernel controls how quickly disorder settles in. 


(22) 
(23) 


Fig. 1. We can use weight-space sampling to directly estimate the naive 
form of the ensemble average light transport in Eq. (8) for stationary kernels. 
We first sample a collection of basis functions based on the kernel and then 
a realization by choosing weights for each basis function (left). This gives us 
a fixed global realization that we can evaluate at arbitrary points in constant 
time. We then use affine arithmetic-based ray tracing [Sharp and Jacobson 
2022] to find ray-surface intersection (right). 


Thin-plate. Williams and Rasmussen [2006] suggest the “thin- 
plate” covariance. It is motivated by the fact that posterior GPISes 
should not “return to zero” for locations far away from the condi- 
tioning points. Instead, it is a more reasonable assumption that the 
magnitude of the sampled realizations should increase as you move 
away from the surface sample. 


kTP (x,y) = 0° (2||x — yll? + 3R||x — yll? +R?) (24) 
302 
Pios 45 (25) 


Note that this covariance is only p.s.d. for ||x — y|| < R, so R has to 
be chosen to be at least as large as the size of the domain. 


3 ALGORITHMS 


In this section, we discuss the rendering algorithms we use to vi- 
sualize Gaussian process implicit surfaces in more detail. These 
algorithms are, essentially, direct Monte Carlo estimators of Eq. (8) 
and Eq. (14). In particular, we formulate them in such a way that we 
avoid ever having to evaluate the GPIS T(t, n) density directly. 


3.1 Global realization sampling via weight-space GPs 


We first describe a method to directly compute Eq. (9). Recall that 
there, we require realizations that span the whole space we want to 
trace paths through at once, and sampling from high dimensional 
multivariate Gaussian distribution quickly becomes expensive, as 
discussed in Section 3. If we limit ourselves to stationary kernels, we 
can apply the weight space decomposition to the Gaussian process 
instead [Wilson et al. 2020]. Here we approximate realizations as a 
weighted sum of a finite number M of basis functions, e.g. random 
Fourier features ¢;(x) = V?/m cos(07 x + ti) where 0; are chosen 
according to the spectral density S(s) = f k(r)Jaj2-1(27rs) ril? dr 
of the covariance function, with Jg/2—1 being a Bessel function of 
order d/2 — 1. Additionally, r; ~ Uniform(0, 27) and d is the dimen- 
sionality of the domain. Sampling a realization is then equivalent 
to sampling the weights for the basis function, and for an uncondi- 
tioned process, the weights are independently Gaussian distributed. 
The realization can then be evaluated at any set of points X in O(N) 


time as 
M 
F(X) =D) wigi(X), wi ~ N(0,1). 


i=1 


(26) 
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This lets us draw correlated samples at a cost only dependent on 
the chosen number of basis functions instead of on the number of 
correlated samples. 


3.1.1 Practical considerations. In particular, once we have sampled 
weights for the basis functions, we can treat the resulting linear 
combination like any other implicit function and trace against it 
by finding the first zero crossing along each ray. Unfortunately, we 
can’t use a fast root-finding method like sphere tracing because 
the Lipschitz constant of our basis grows linearly with the number 
of basis functions, forcing sphere tracing to take very small steps. 
Instead, we rely on the slower but more general method of root- 
finding via affine arithmetic [Sharp and Jacobson 2022]. Because the 
basis functions use a limited number of operations, implementing 
them in affine arithmetic is relatively straightforward. When we 
have found an intersection, we can analytically compute the nor- 
mal at that point and do shading as for any other implicit surface. 
We illustrate this in Fig. 1. Note that because the basis functions 
use non-linear operations (cosines), the bounds produced by affine 
arithmetic are not tight. Nonetheless, we are guaranteed not to miss 
intersections. Weight-space GPIS rendering makes it practical to 
get results equivalent to the Global © model, assuming that the 
number of basis functions is chosen to be large enough. This means 
that this method is mainly useful for verifying more general ones, 
such as the function-space algorithms we will discuss next, in sim- 
plified scenarios. There are several practical concerns when using 
weight-space GPs for GPIS rendering. The most obvious one is that 
the standard weight-space formulation can only support stationary 
covariance kernels. The mean can still vary arbitrarily, allowing us 
to model interesting scenes, but we can not vary properties such 
as the roughness of surface-type GPIS or include surface-type and 
volume-type GPIS in the same scene. Choosing the correct number 
of basis functions is critical, and the exact number depends on the 
covariance kernel. When choosing a low number of basis functions, 
the underlying kernel is approximated poorly, which leads to sig- 
nificantly different realizations (see Fig. 2). Empirically, we observe 
that kernels with longer-than-exponential tails, such as the rational 
quadratic kernel, require a much larger number of basis functions 


x 


| . function-space +- 


Fig. 2. When we choose the number of basis functions we need to take the 
roughness of the covariance kernel into account (top: smooth SE kernel, 
bottom: rough RQ kernel). For rougher processes, too few basis functions 
result in biased realization samples that do not predict the ground truth 
kernel well. (left to right, M = 10, 100, 1000, ground truth). 


than, for example, the squared exponential kernel. Finally, one has to 
be able to sample from the spectral density of the chosen covariance 
kernel. This is not always trivial to do, and we provide procedures 
for the stationary kernels used in this work in Section 2.1 of the 
supplemental. 


3.2 Practical function-space sampling strategies 


In Section 4.2 of the main text we discussed our function-space 
sampling strategy and left off with two practical issues to solve 
in an implementation: Determining scene bounds along rays and 
limiting the number of correlated values that have to be sampled 
jointly. Here we go into more detail on how we handle both of these 
issues. 


Determining scene bounds. One simple and practical solution to 
define scene bounds is to prescribe them a priori. For example, we 
can simply intersect the ray with a box or a sphere that we define 
as the limits of the scene, resulting in only having to sample a 
realization over a finite ray segment. A more principled option is to 
consider that many scenes are naturally finite in extent. In our case, 
this occurs when u(x) — oo for large x, e.g. when the mean surface 
is modeled using an SDF. Then, for large x, P(f (x) < 0) — 0, i.e. the 
probability of sampling a value that we see as “inside” a GPIS goes 
towards 0. We can define the bounds of the scene as the region for 
which P( f(x) < 0) > e for some small epsilon. That is, we discard 
regions of space for which we are very unlikely to sample a zero 
crossing. Evaluating this “point-wise occupancy” is cheap as it is 
just the CDF of a mono-variate normal distribution. We can intersect 
the ray with the level set {x | P(f (x) < 0) = e} using ray marching. 
Accuracy here is not critical as long as we are conservative. This 
gives us a finite ray segment to work with as long as our scene is 
finite. The only case this does not cover is if we have an infinite 
scene. The most common example of this would be a scene filled 
with an infinite volume-type GPIS. Here we will have to rely fully 
on the second level of our strategy. 


Progressive ray segment sampling. There are two cases where we 
can’t simply distribute sample points along a ray segment: We have 
an infinitely large scene, or the ray segment is long compared to the 
length scale of the covariance kernel (in which case the O(n?) scal- 
ing of function-space sampling would make sampling a sufficiently 
dense set of points too slow). In practice, we found that taking 
anything more than 256 correlated samples at a time significantly 
impacts render times. To overcome this, we apply the same strategy 
we applied in Section 4 to path segments to ray segments. That is, 
we sample realizations Sx1;%t41) ~ GP (x6, Xti1) ICAL! with tọ = 0 
and tj41 = ti +n- At, where At is a step size chosen based on the 
covariance kernel, and ¢’ the condition based on the chosen mem- 


ory model. We repeat this until I (ey, = 0, ie. Fx; tit) 
contains a zero crossing. We then proceed as before, sampling a 


normal followed by an outgoing direction. This process is illustrated 
in Fig. 8. 
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3.3. Thoughts on next-event estimation 


Our function-space sampling approach enables us to employ next- 
event estimation to a much greater degree than the global weight- 
space sampling approach. When determining a global realization 
before tracing a path, we fix not only the realization’s values but 
also its gradients and, hence, normals. Then, during tracing, there is 
exactly one possible normal at each intersection point. This means 
that, when using a mirror micro-BSDF, we cannot perform next- 
event estimation since, with the normal at the intersection point 
being deterministic, the outgoing ray direction is completely deter- 
mined. In function-space sampling, the normal at an intersection 
is not determined ahead of time. Instead, we sample it from a dis- 
tribution conditioned on the realization values we saw along the 
incoming ray. This realization only fixes the directional derivative of 
the GP at the intersection point but leaves us two additional degrees 
of freedom to choose a gradient and, hence, a normal. Assuming 
non-degenerate covariance kernels (i.e. volume-type GPISes, that 
do not assume perfect correlation along any axis), the resulting 
gradient distribution is non-degenerate and assigns some density to 
the whole hemisphere of normals facing the ray. Hence, even when 
using a mirror micro-BSDF, we can always find a normal such that 
the reflected ray points into the desired NEE direction. Computing 
the density of this distribution for a given normal is not trivial. The 
conditioned distribution of gradients is a 3D multivariate Gaussian 
distribution (no matter how complex the conditioning is). To com- 
pute the (unnormalized) density of sampling a given normal, we 
simply need to integrate over all possible gradients that normalize 
to that normal as 


p(n) = f fani (27) 


where f is the pdf of the conditioned gradient distribution at the 
intersection point. Our preliminary investigation has shown that 
the integral in Eq. (27) is tractable analytically, and we should be 
able to compute the normalization constant to turn p(n) into a true 
pdf. This would then allow us to evaluate the conditioned visible 
normal distribution, enabling (single scattering) NEE for mirror 
micro-BSDFs. We leave this for future work but still apply classical 
NEE after sampling a normal when the micro-BSDF is diffuse. 


3.4 Implementation 


We implemented the sampling strategies described in this section 
in the Tungsten renderer [Bitterli 2018] with a focus on correctness 
over performance. A GPIS is treated as a participating medium that 
allows for sampling free-flight distances and computing transmit- 
tance, and provides a phase function at scattering locations. Light 
and camera positions are uncorrelated from the GPIS and we often 
place uncorrelated, non-GPIS surfaces in the scene alongside the 
GPIS we are investigating. We could, of course, represent these as 
uncorrelated GPIS as well, but that would unnecessarily increase 
rendering times without aiding in the validation or understanding 
of our method. We use standard data structures such as OpenVDB 
grids to store volumetric scene data (variance, length scale, and 
mean). For the mean, we also support using a mesh directly and 
computing SDFs on the fly. 


4 UNCERTAINTY QUANTIFICATION USING DOUBLE 
MONTE CARLO SAMPLING 


We have a model M(x) with uncertain parameters ® ~ po. We 
would like to compute the moments of the distribution of model 
outputs qx(yx = M?(x)). This problem is known as uncertainty 
quantification. If we have access to a deterministic method of evalu- 
ation M® (x), we can compute 


Eo~po [M°(x)| = J M? (x) dpa (®) (28) 
Vo-pe |M] = Eo~po [M*(2)?] -Eo-p [MPO] 9) 


2 
i MP? apace) -| f M(x) apa0) 
(30) 


and use Monte Carlo integration and sample variance to compute 
an estimate of the mean and variance, respectively as 


Eo~py [M®(x)] y 5 M* (x) Di~ po (31) 
—— ee i 

Vonpe [MP0] = 3 DMP)? (DMM)? di ~ Po 
i=0 i=0 

(32) 


Note that due to the square in the second term and Jensen’s inequal- 
ity, we need a relatively large number of samples to get an unbiased 
estimate of variance. In practice, this is not a big issue since we 
tend to have the samples available anyway if we want to compute 
a relatively converged estimate of the mean. That is, if we have 
enough samples to estimate the mean, we tend also to have enough 
samples to estimate variance. 

Unfortunately, in many graphics applications like rendering, the 
model itself is often in the form of a complex integral equation, such 
as 


M’) = J m? (x, y) dy (33) 


and has to be estimated using Monte Carlo techniques. That is, we 
only have access to an estimate 


TIN S m? (x yi) 
M(x) =) oi ~ Py (34) 
2 pyi) °F 
This is not an issue when computing the mean of our model predic- 
tions. Since a one-sample Monte Carlo estimator is unbiased, that 
is 


ee ® 

M?(x) =E [M] =Ey-p, [2], (35) 
we can write 
Eo~pa [M°(x)| = Eo~py [E [MPC 1e] (36) 
P(x,y) 
sana pen eR] on 
Cy) 

m Ei ar | (8) 
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and then use Monte Carlo integration to estimate the expectation as 


m® (x. y) 1 mt: (x, yi) 

E~ = . N : ®; ~ pa, Yi ~ 
manne o| pyly) nd ast, PSHP 

(39) 


Note that even though we are sampling two random variables now, 

we can still just average N independent evaluations of m® (x, yj). 

This is one of the central benefits of Monte Carlo integration; its 

convergence does not depend on the dimensionality of the problem. 
But for the variance, we have 


m® (x, y) 


Py) i) 


Vo~pe [MPoo| # Vo~po.y~py | 
Intuitively, this is because V [MG N| ol + 0 for N < œ and we 


have, according to the law of total variance 


von [iy] - 
Eyi, YN~Py [Vo~pe [MPC | yo „un | 
+ Vy s yN~py [Bo~ps [MPC lyst un | . (41) 


Here, we can see that for finite N, the total variance is comprised 
of the expected value of the variance of our model due to model 
parameters and the variance due to the Monte Carlo estimation of 
our model. 

We only recover the ground variance of the model due to the 
model parameters as N — ov: 


Nim Ey yn~py [Vo-po [MOn Iyo: 9] = Voxpo [M] 
(42) 
Aim Vy yn~py [Bo~po [MPG Lynsay |] = 3) 


Intuitively, now that means that if we want to compute the sample 
variance of our model due to the model parameters unaffected by 
variance due to Monte Carlo estimation of the model, we need to 
sample many yj ~ py for each 0; ~ po. Even if we reuse the same 


set of ys, we still need to evaluate M?! (x) y separately for each 0; 
and compute sample variance as 


N K 


oa 1 1 m?i (x, yj) 
Vo~pe MP) | nk = N dK >, U y 
i= j= 


1 aa mi (x, yj) 2 
-9 Gi ~ poyj~ py (44) 


This now results in quadratic complexity O(NK) where N controls 
the bias due to the limited number of parameter samples and K 
controls the bias due to the additional variance in the Monte Carlo 
estimator. Hence, while it is certainly possible to do uncertainty 
quantification via Monte Carlo rendering of GPISes, it is not trivial 
to do so efficiently and we leave this for future work. 


REFERENCES 


Benedikt Bitterli. 2018. Tungsten Renderer. https://github.com/tunabrain/tungsten/ 

Christian Grosche and Frank Steiner. 1998. Handbook of Feynman path integrals (1998 
ed.). Springer, Berlin, Germany. https://doi.org/10.1007/bfb0109520 

Nicholas Sharp and Alec Jacobson. 2022. Spelunking the Deep: Guaranteed Queries on 
General Neural Implicit Surfaces via Range Analysis. ACM Transactions on Graphics 
(Proceedings of SIGGRAPH) 41, 4 (July 2022), 107:1-107:16. https://doi.org/10/grnsz3 

Christopher KI Williams and Carl Edward Rasmussen. 2006. Gaussian processes for 
machine learning. Vol. 2. MIT press Cambridge, MA. 

James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc 
Deisenroth. 2020. Efficiently Sampling Functions from Gaussian Process Posteriors. 
In Proceedings of the 37th International Conference on Machine Learning. PMLR, 
10292-10302. https://proceedings.mlr.press/v119/wilson20a.html 


From microfacets to participating media: A unified theory of light transport with stochastic geometry — Supplemental + 7 


kernel & derivatives f.s. realization w.s. realization posterior occupancy 


(oa 
fa 7 F P not implemented 
re 


ga 
as 
we 
—2.5 
Q 
ee 
Se 
—2.5 
A 
n 
2 
Ou 4 
Fue not implemented 


Fig. 3. We show an overview of some of the kernels that we implemented for our method. To give an intuition for the family of GPISes each kernel produces, 
we show samples from the prior (produced via both function-space and, where applicable, weight-space methods) and posterior occupancy. 


